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Abstract 



In the field of numerical integration, methods specially tuned on os- 
cillating functions, are of great practical importance. Such methods are 
^. ■ needed in various branches of natural sciences, particularly in physics, 

^^ ' since a lot of physical phenomena exhibit a pronounced oscillatory behav- 

ff") , ior. Among others, probably the most important tool used to construct 

Q\ • efficient methods for oscillatory problems is the exponential (trigonomet- 

Cn ' ric) fitting. The basic characteristic of these methods is that their phase 

r — , lag vanishes at a predefined frequency. In this work, we introduce a new 

f*^ ■ tool which improves the behavior of exponentially fitted numerical meth- 

QiQ | ods. The new technique is based on the vanishing of the first derivatives of 

^^ . the phase lag function at the fitted frequency. It is proved in the text that 

these methods present improved characteristics in oscillatory problems. 



PACS: 0.260, 95.10.E 

1 Introduction 

The numerical integration of systems of ordinary differential equations with 
oscillatory solutions has been the subject of research during the past decades. 
This type of ODEs is often met in real problems, like the Schrodinger equation 
and the N-body problem. For problems having highly oscillatory solutions, the 
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standard non-specialized methods can require a huge number of steps to track 
the oscillations. One way to obtain a more efficient integration process is to 
construct numerical methods with an increased algebraic order, although the 
implementation of high algebraic order methods is not evident. 

On the other hand, there are some special techniques for optimizing nu- 
merical methods. Trigonometrical fitting and phase-fitting are some of them, 
producing methods with variable coefficients, which depend onu = uih, where w 
is the dominant frequency of the problem and h is the step length of integration. 
More precisely, the coefficients of a general linear method are found from the 
requirement that it integrates exactly powers up to degree p+1. For problems 
having oscillatory solutions, more efficient methods are obtained when they are 
exact for every linear combination of functions from the reference set 

{l,x,...,x K ,e ± ^,...,x p e ± ^} (1) 

This technique is known as exponential (or trigonometric if \x ~ iu) fitting and 
has a long history [3], [TJ]. The set ([T]) is characterized by two integer param- 
eters, K and P . The set in which there is no classical component is identified 
by K = —1 while the set in which there is no exponential fitting component 
(the classical case) is identified by P = — 1. Parameter P will be called the 
level of tuning. An important property of exponential fitted algorithms is that 
they tend to the classical ones when the involved frequencies tend to zero, a fact 
which allows to say that exponential fitting represents a natural extension of the 
classical polynomial fitting. The examination of the convergence of exponential 
fitted multistep methods is included in Lyche's theory [14]. There is a large 
number of significant methods presented with high practical importance thats 
have been presented in the bibliography (see for example [23], [8], [17], pQ, [2], 
0], [S], HE], 0, HD], 02], 13 , 02], [22], [2S]- The general theory is presented in 
detail in [TT] . 

Considering the accuracy of a method when solving oscillatory problems, it 
is more appropriate to work with the phase-lag, rather than the usually used 
principal local truncation error. We mention the pioneering paper of Brusa 
and Nigro [5J, in which the phase-lag property was introduced. This is actually 
another type of a truncation error, i.e. the angle between the analytical solution 
and the numerical solution. On the other hand, exponential fitting is accurate 
only when a good estimate of the dominant frequency of the solution is known 
in advance. This means that in practice, if a small change in the dominant 
frequency is introduced, the efficiency of the method can be dramatically altered. 
It is well known, that for equations similar to the harmonic oscillator, the most 
efficient exponentially fitted methods are those with the highest tuning level. 
In the case of the Schrodinger equation, this result was already obtained for 
particular two- and four-step exponentially fitted multistep methods based on 
an expensive error analysis, see for example [12], [IE], [32] and [20] . 

In this paper we present a methodology for optimizing numerical methods, 
through the use of phase-lag function and its derivatives with respect to v. 
More specifically, given a classical (that is with constant coefficients) numerical 



method, we can provide a family of optimized methods, each of which has zero 
phase lag (the case of trigonometric fitting) or zero PL and PL' or zero PL, 
PL' and PL" etc. With this new technique we provide methods that perform 
well during the integration of the Schrodinger equation for high values of energy, 
but also that perform well on other real problems with oscillatory solution, like 
the N-body problem. 

2 Phase lag analysis 

Below we will consider for simplicity only first order differential equations, al- 
though the same results can be easily obtained for second order equations too. 
Consider the test problem 

' l! ' [h) =iu; y(t), 2/(0) = 1 (2) 



dt 



with exact solution 



y{t) = e"-«* (3) 



where wo is a non-negative real value. Let 3>{h) be a numerical map which when 
it is applied to a set of known past values, it produces a numerical estimation of 
y(t + h). If we assume that all past values are known exactly, then the numerical 
estimation y(t + h) of y(t + h) will be 

y(t + h) = a{uj h) ■ e ^ot+H^oh)) ^ 

while the exact solutiion is e i ("°*+ w °' 1 ). Then 

L = y^lR = a(w h)e- i ^ h -^° h ^ (5) 

y(t + h) 

In the above equation l|5]). the term a(u>oh) is called the amplification factor, 
while the term 1(u)qK) = uj^h — <j>(u>oh) is called the phase lag of the numerical 
map. In the case that a(uj Q h) = 1 and I(luqIi) = 0, we say that the numerical 
map $(/i) is exponentially fitted at the frequency loq and at the step size h. 

Suppose now that the method $(/i) has been designed in order to solve 
exactly equation ^. But in practice, only an estimation of the frequency u>o is 
known. Thus, it is of great importance to know the behavior of the method at 
frequencies close to the estimated one, so we apply the method to the equation 

^ = &*(*), y(0) = 1 (6) 

and calculate the phase lag l(u), where u — uh. Since the method integrates 
exactly equation (|2|), the phase lag function l(u) has a zero at point uq = u)oh 
and is given by 

l(u)=u-(f>(u) (7) 



But since 4>(uq) — uq we have 

l{u) =u-(f>{uo) ^~i»="o - u o) + 2^ ^ M » l"="o rj ( 8 ) 

or 

„ , , J dd>(u). \ ^d^Mu), («-«o) n 

v ' n— 2 

and thus we conclude that in order to maximize the phase lag order at frequen- 
cies close to mo = ojoh we must at least have 

^1 - 1 (101 

while in order to obtain even higher order in phase lag, the higher derivatives 
of <f>(u) at point uq must vanish. Since 

l'(u) = 1 - <t>'(u), l {p) (u) = (t> {p) (u), p > 1 (11) 

we have the following theorem 

Theorem 1 . Consider a linear method <t which solves exactly the equation 
(|2j) and when it is applied to the equation y' = iuy, with u) = ujq + 8, it produces 
a phase lag function l(u), u = Luh. Then, if the phase lag function has its s first 
derivatives at point vq = uj^h equal to zero, then the phase lag function l(u) is 
of order at least s in 8. 

Supose now that the method $ depends on M independent parameters and 
let $ c the classical method which is constructed by setting M equations max- 
imizing the algebraic order of the method. Let $ s be the method which is 
constructed with M — 1 — s equations maximizing the algebraic order, 1 equa- 
tion for vanishing the phase lag at a frequency ui^h and s equations for vanishing 
the s first derivatives of the phase lag at the same point uj^h. It is easily now 
calculated that the local truncation error of the method <£> s is given by 



lie* = 



«(u) (e- a{u) - l) (12) 



when the working frequency ujq — ► 0. But since the method is at least of order 
M — 1 — s we have 

lte a = h M - s (l(u )-il'(uo)(w-u )h 

+ Er= 2 (A J ^ ) M + g J a (1) M,...,^- 1) M))^' UdJ 

with gj be a polynomial function of the j — 1 first derivatives of l(u) at point uq 
with <7j(0, 0, . . . , 0) = 0. Since now l(u ) = and the s first derivatives of l(u) 
at uo vanish, we have that in the limit luq — > 

lte s = ch M+1 = lte c (14) 

where lte c is the local truncation error of the classical method $ c . Thus we 
have the following theorem 



Theorem 2. Consider a linear method $ s which is constructed by de- 
manding (i) maximal algebraic order and (ii) that the phase lag and its s first 
derivatives vanish at some given frequency wq. Then, when ojq — > 0, the method 
is identical with the classical one & c which is constructed by demanding only 
maximal algebraic order. 

3 Numerical results 

In order to follow the dynamics of the method constructed, vanishing the first 
derivatives of the phase lag function, we consider the simple 2 — step symmetric 
formula 

y n _i - 2 • y„ + y n+1 = h 2 (b ■ (f n -i + f n +t) + hfn) (15) 

for the solution of the 2nd — order equation 

%-m tie) 

The coefficients of the method are calculated in three different cases as follows: 

1. b c are the coefficients for the classical method, where only the maximiza- 
tion of the algebraic order is taken into account 

2. b* are the coefficients for the method where trigonometric fitting in a 
frequency u> (y = u>h) is taken into account and 

3. 6 s are the coefficients for the method where both trigonometric fitting 
and vanishing of the first derivative of the phase lag function is taken into 
account. 

4. b sd are the coefficients for the method where both trigonometric fitting 
and vanishing of the first and second derivatives of the phase lag function 
is taken into account. In order to obtain this, the coefficient of y n in 
equation (TT5]is perturbed from —2 to —2 + a(v). 

The coefficients for the three cases are given 

K - 1, K -| d7) 

,t 1 1 2 17 a 31 fi 691 o 

bi = — H v 2 + v 4 H v 6 -\ v 8 + 

12 120 20160 362880 79833600 

5461 V" + 929569 „» + Q(i>") (18) 



6227020800 10461394944000 
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-v 
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3113510400 5230697472000 



« 12 + (v 14 ) (19) 



bo 
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usd 1 1 2 41 4 121 9 6 8887 8 
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+ 0(v 13 ) (22) 



&f 



5 1 , 17 4 1811 fi 13817 o 

6 40 2016 1814400 79833600 

12478951 10 24838031 12 , ir . 

v 10 H u 12 + 0(u 13 ) (23) 

435891456000 5230697472000 v ; v ' 



a(v) = -—v 6 - —v s —v w ^^12 _ 

v ; 240 2016 11520 159667200 

62879 v" + Q(„") (24) 

26417664000 v ; v ' 

The problem under test is the 2-body problem and the 2nd order equation 



is 

d 2 ( Vi{t) \ _ 1 ( Vi(t) 



dt 2 \ y 2 (t) J ~ y/yl(t) + y*(t) V 2/2(4) 
with initial conditions 

I&(0) \_( l-e\ ( yi (Q) 
2/ 2 (0) J I J ' I y 2 (0) 



and exact solution 



and 



2/i (*) ^\ / cos(it) -e 



y 2 (i) ; - 1 Vi 



e 2 sin(u) 



(25) 



(26) 



(27) 



u - €sin(u) - t = (28) 

In Figure Q] the error in position calculation is plotted for the four methods. 
The working frequency in the trigonometric fitted methods has been estimated 
by u> = — j- [I] . Finally in Figure [2] the phase lag for the four methods is 

shown as a function of the frequency. 



4 Conclusions 

A new technique has been developed in order to improve the behavior of ex- 
ponential (trigonometric) fitted numerical methods for the integration of os- 
cillatory problems. The new technique is based on the vanishing of the first 
derivatives of the phase lag function, thus decreasing the sensitivity of the nu- 
merical method to frequency variations. Moreover, it has been shown that the 
new method becomes the classical one (the one who maximizes the algebraic 
order) when the working frequency tends to zero. 
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List of Figures 

1. Mean error in position calculation in the two body problem for eccentricity 
e = 0.5 and step size h = 0.1. (o) is the algebraic fitted method, (H) is 
the trigonometric fitted method, (x) is the trigonometric fitted method 
with the first derivative of the phase lag function equal to zero and (<)) is 
the trigonometric fitted method with both first and second derivatives of 
the phase lag function equal to zero. 

2. The phase lag of the four methods as a function of frequency (o) is the 
algebraic fitted method, (□) is the trigonometric fitted method, (x) is 
the trigonometric fitted method with the first derivative of the phase lag 
function equal to zero and (<0>) is the trigonometric fitted method with 
both first and second derivatives of the phase lag function equal to zero. 
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Figure 1: Figure 1. 




0.04 0.06 0.08 0.1 0.12 0.14 

frequency (x time step) 



0.16 



Figure 2: Figure 2. 
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